Solving for three-dimensional central potentials using matrix mechanics 
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Matrix mechanics is an important component of an undergraduate education in quantum me- 
chanics. In this paper we present several examples of the use of matrix mechanics to solve for a 
number of three dimensional problems involving central forces. These include examples with which 
the student is familiar, such as the Coulomb interaction. In this case we obtain excellent agreement 
with exact analytical methods. More importantly, other interesting 'non-solvable' examples, such 
as the Yukawa potential, can be solved as well. Much less mathematical expertise is required for 
these methods, while some minimal familiarity with the usage of numerical diagonalization software 
is necessary. 



I. INTRODUCTION 

An important component of the undergraduate train- 
ing in quantum mechanics is the solution of the three 
dimensional Schrodinger equation for a particle that ex- 
periences a central potential, i.e. one in which the po- 
tential is a function only of the distance from the origin. 
Then, the Schrodinger equation is separable in spherical 
coordinates, and the angular part is determined analyt- 
ically, in terms of the well known spherical harmonics 
and the associated Legendre polynomials. 1 What then 
remains is the solution of the radial part of the wave 
function, and typically some examples are worked out, 
like the infinite spherical well, the 3D harmonic oscilla- 
tor, and of course, the hydrogen atom. 

The solution to the radial part of the wave function 
usually requires somewhat advanced mathematics, and 
tends to go in one of two ways, either (i) solution by 
recognition, or (ii) by power series. The first method 
amounts to declaring that the radial differential equa- 
tion is one that has been studied for more than a century, 
and is 'easily' recognizable as the 'insert famous name' 
Equation, and therefore has 'insert famous name' func- 
tions as solutions. Even worse from a student's point of 
view, is to declare that the solution is a confluent hy- 
pergeometric function with appropriate arguments, and 
perhaps leave as an exercise which famous name is asso- 
ciated with which arguments. The second method is a 
little more satisfying, in that the student ends up con- 
structing the function that turns out to have a famous 
name associated with it, but there are often a number of 
preliminary steps required, whereby the asymptotic be- 
haviours are 'peeled off', so that what remains is a simple 
polynomial; this procedure is straightforward to those 
that are familiar with it, but to a novice in both quan- 
tum mechanics and in differential equations, the process 
can be somewhat daunting. 

This level of mathematics has its benefits, and is of 
course a required component of a physicist's toolkit. 
However, at this stage of a student's career it can also 
serve to dampen their enthusiasm for physics. As educa- 
tors it is also important to consider that for students who 
eventually do not pursue a career in physics (or math- 
ematics), extensive knowledge of Laguerre polynomials 



will probably not help them in their future career. Fur- 
thermore, this way of solving problems will not be too 
helpful when it comes to examining potentials that are 
not tractable by either of these methods. 

At the same time, matrix mechanics is generally 
taught only in the abstract, with real implementations 
relegated to more advanced degrees, and usually in the 
context of many-body physics. We therefore suggest a 
general purpose numerical method for solving problems 
involving central forces in three dimensions; the results 
obtained are necessarily approximate but very accurate. 
While this method should not replace the teaching of 
exact analytical methods referenced above, it provides a 
tool for learning the use of matrix mechanics methods, 
and for understanding the behaviour of the solutions of 
various potentials that cannot be solved analytically. It 
has even been used very recently to provide insight to 
the differences between the 2s and 2p electrons in atomic 
lithium. 2 This method follows that already developed for 
one dimensional potentials in Ref. 3. The virtue of this 
approach is that it requires only linear algebra and in- 
tegral calculus, topics normally covered in a student's 
first year of university studies. The difficult part is that 
students need to have access to software tools at some 
basic level to carry out the linear algebra, and, in some 
cases, to perform the integrations that are required. Our 
experience has been that this part is difficult for some 
students; however, it is our belief that some familiar- 
ity with Matlab, or Maple, or Mathematica will have 
broader application for the average student in the long 
run than a knowledge of non-elementary functions. 

We begin with an example which can be first addressed 
by standard methods, the Coulomb potential, specifi- 
cally for the hydrogen atom. In this way students can 
readily check their answers. The Coulomb potential hap- 
pens to be one of the most difficult examples to use, 
however. Because of its long range it supports an infi- 
nite number of bound states; as will be shown below it 
is impossible to recover all of these through the present 
approach, but a careful study of this problem will help to 
highlight the limitations and subtleties of this approach. 
Probably a few tricks could be adopted to circumvent 
this difficulty, but this would run counter to our goals, 
as this method should be generally applicable to any po- 



tential that supports bound states. 

We will then examine the finite spherical well, and 
determine for example, the critical depth required for 
at least one bound state to exist. This is also known 
analytically, and so will provide a benchmark for the 
present method. 

Finally, we will examine solutions for the so-called 
Yukawa potential, a useful potential both in nuclear 
physics, where it was used to model meson exchange be- 
tween nucleons, and in condensed matter physics, where 
it is used to model Coulomb interactions whose range 
has been shortened due to screening. Solutions for this 
potential generally require advanced applications of per- 
turbative or variational methods. We will make compar- 
isons of our results with these. 

We should emphasize that we will not comment on 
numerical methods in this paper. We assume that stu- 
dents have access to software that can compute desired 
integrals and diagonalize reasonably large matrices. It is 
assumed that the necessary training for this procedure 
is a prerequisite to a student's first course in quantum 
mechanics^ 



II. FORMULATION OF THE MATRIX 
MECHANICS PROBLEM 

For a central potential the solution to the time- 
independent Schrodinger equation is separable, 



iP(r,e,<p)=R(r)Y t m (e, 



(1) 



where we have already written down the solution to the 
angular part — it consists of the spherical harmonics, 
which are functions of the standard spherical angles. Fol- 
lowing the usual procedure, one can replace the radial 
function R(r) with u(r) = rR(r), and arrive at the so- 
called radial equation, 



h 2 d 2 u{r) 
2mn dr 2 



+ V eS (r)u(r) = Eu{r) 



(2) 



which is identical to the one dimensional Schrodinger 
equation for a particle of mass mo, except that the effec- 
tive potential contains an additional so-called centrifugal 
term, 



Veff(r) = V{r) 



H 2 £(£+1) 
2mn r 2 



(3) 



Note also that J °° dr\u(r)\ 2 = 1, and that r ranges from 
to oo. In addition, because this Schrodinger equation 
is for u(r) = rR(r), and R(r) is well-behaved at the 
origin, then a boundary condition is that u(r = 0) = 0. 5 
Following Ref. we embed this potential in an infinite 
square well^ extending from r = to r = a, where a is 
some cutoff radius, whose value will influence the results 
in a manner to be explained below. 

Figure 1 shows the examples of the potentials we will 
use in this paper, plotted along with the infinite square 



■Coulomb + Infinite Well 
Yukawa + Infinite Well 
Spherical Well+ Infinite Well 

■Arb. Potential + Infinite Well 




FIG. 1. (color online) A plot of the various potentials to 
be used in this paper, along with the infinite square well 
"embedding" potential placed between r = and r — a. Note 
that we have used ao/a — 1/10 in the figure, and b = a/10 
for the finite spherical well. Also shown is an arbitrarily 
complicated potential well, to illustrate the point that this 
potential poses no further difficulty compared to the others, 
with this method. 



well "embedding" potential. We also include an arbi- 
trarily complicated potential well shape, to emphasize 
that this method can solve for any such potential. The 
rationale for this choice is that the embedding potential 
allows for a simple set of basis states, which are simply 
the eigenstates of the infinite potential well, 



<l>n(r) 



2 , mrr 
sm ), 

si 



with eigenvalues 



E° = 



n 2 h 2 n 2 
2rriQa 2 



(4) 



(5) 



The embedding potential enforces that the function is 
zero at the origin, u(r = 0) = 0, but also now requires 
the function to vanish at the other wall, u(r = a) = 0. 
Such a formulation will work reasonably well for bound 
states in attractive potentials, and provides a basis set 
that is most familiar to students. Now we write the 
radial equation in Dirac notation as 



[H + V cS ]\u)=E\u) 



(6) 



where Ho includes both the kinetic energy and the infi- 
nite square well potential, so that 



Hi 



0\9n 



(7) 
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If we now expand \u) in terms of this basis set, 



oc 



i\<t>m), 



(8) 



and substitute this into Eq. (j6|), followed by an inner 
product with each bra (<J) n | , we obtain the matrix equa- 
tion 



oo 

E 



(9) 



where the matrix elements are given by 



(<j> n \(H + V eS )\<j>„ 



S E v 



dr sin 



, nirr 



){V(r) 



h 2 £(£+l) 
2mor 2 



} sin ( — )(ao) 



and S nm is the Kronecker delta function. 

Eq. (fTUf is readily evaluated for any bound state po- 
tential, numerically if need be. We will begin with the 
Coulomb potential experienced by an electron of reduced 
mass toq near a positively charged proton, 



Vboul(r) 



47reo r 



(11) 



where =Fe is the charge of the electron (proton) and eo 
is the vacuum dielectric constant. Note that both the 
Coulomb potential and the centrifugal term are singular 
at the origin, but that the integrand in Eq. (fTU|) is not, 
and varies smoothly, apart from the oscillations of the 
sine functions. 

Two issues should be raised before we proceed; one is 
that the matrix size in the eigenvalue problem posed in 
Eq. ([9]) is infinite. This will be dealt with by utilizing an 
upper cutoff n max , and increasing the value of this cut- 
off until the results are converged. For large quantum 
numbers, the basis states exhibit two related properties; 
first, they have increasing energy, and for this reason, 
may be viewed as more and more irrelevant for contri- 
butions to low energy states. However, concomitantly 
they have finer spatial resolution. Often it is the finer 
spatial resolution that is required to accurately describe 
a low-lying state, so sometimes a larger basis is required 
than one might think if only energy considerations are 
used. Either way, numerical convergence is attained once 
basis states with very high energy and very sharp spa- 
tial resolution are not needed to describe the problem at 
hand. 

The second issue concerns the value of the width of the 
well, a, or equivalently, the cutoff in radial distance. The 
natural unit of distance in the Coulomb problem is the 
Bohr radius, ag = (Anen/ e 2 ){h 2 /to), since we know in 
advance that for the Coulomb problem the bound states 
decay exponentially with radial distance r. As we shall 
see, however, exponential decay is not as strong as we 
would like, particularly when there is a numerical coeffi- 
cient in the exponent that stretches out the exponential 
decay, so that a cutoff in the radial distance is required 



to be many times (10-20) the Bohr radius to get very 
accurate results for the ground state. For excited states 
this cutoff will have to be higher, to achieve the same 
level of accuracy. 

For large well widths, however, the basis states dete- 
riorate in spatial resolution. To achieve the same spatial 
resolution, therefore, we will need to increase the value of 
n m ax- We will return to these comments as we examine 
particular examples in the following sections. 



III. THE COULOMB POTENTIAL 

For the Coulomb potential we will focus on I = 0, to 
illustrate the method. Note that the integral in Eq. (flO|) 
can be written in terms of the so-called cosine integral. 
However, in the spirit of avoiding non-elementary func- 
tions (this one is more easily evaluated in the form of 
the integral written in Eq. (jlOl) anyways), we simply do 
the integral numerically. Note that in the interest of cal- 
culating as few of these integrals as is needed ahead of 
time, it is best to use a trigonometric identity^ before 
evaluating the integral. We obtain 



Coui; 



dr sin ( ) — sin ( 

47T£o a Jo a r a 



-2—E Q { L 1 (n + m) -L x (n- 



with 



Lx(m) 



dx 



1 — cos (rmrx) 



. (12) 



(13) 



where we used x = r/a, added and subtracted unity 
to the cosines, and, in the second line of Eq. (fl"2j) . 
we adopted the natural energy unit in the problem, 
E = h 2 /(2m al), which is one Rydberg (« 13.606 eV). 
Similar simplification is applicable for the centrifugal 
term if needed. 

Having decided on a value of a/ do, and a particular 
maximum size for the matrix, n max , it is now a simple 
matter of evaluating the matrix elements and substitut- 
ing into Eq. . We rewrite this equation in dimension- 
less units: 



(14) 



where e = E/Eq, and 



En 



= S nm (^na /ay 



— 2 — \Li(n + m)—Li(n — m)} 
a L ' 

+£(£ + l)(2f) 2 {L 2 (n + to) - L 2 {n - to)}, (15) 



where 



_ , . f , 1 — cos (to7to;) 
L 2 [m) = / dx . 

Jo x 



(16) 
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Figure 2 shows an example of the calculation with 
a /clq = 50 and n max — 200. Several characteristics 
should be noted. The exact results, e eX act = E c ^ act /EQ = 
— 1/n 2 , are given by the square symbols. They are of 
course all negative, though they become more difficult 
to distinguish from zero as n increases. The cross-hairs 



TABLE I. Table 1. Results for a/a = 10. 



T^max 


Ei/E 


E 2 /E 


50 


-0.99915 


-0.22546 


100 


-0.99983 


-0.22558 


200 


-0.99992 


-0.22560 


400 


-0.99993 


-0.22560 



o 
LU 

— I 



1.5 



1.0 



0.5 



0.0 



-0.5 



-1.0 




□ □□□□□□□□□□□□ 



□ -1/n (exact result) 

+ numerical results (with a/a =50 and n nmax = 200) 



10 

n 



15 



20 



FIG. 2. (color online) Energy levels for the Coulomb poten- 
tial vs. the quantum number n (£ = 0). Exact analytical 
results are shown with the squares. Numerical results ob- 
tained with an embedding infinite square well with width 
a — 50ao, where ao is the Bohr radius, are shown with cross- 
hairs. The numerical results reproduce very accurately the 
first four bound state energies. Eventually, the energies be- 
come positive (unbound) due to the embedding potential, and 
for very large quantum number n, they will vary as n 2 (shown 
with a blue dotted curve), suitable for an infinite square well 
of width a. A more accurate result (Eq. (|17[) ~) is shown with a 
red solid curve, and can be derived from perturbation theory. 



represent the result of our numerical calculation. Only 
the first six values of the energy are negative, and of 
these, the first four are very accurate; the remainder be- 
come positive and in fact approach the results expected 
from an infinite square well potential of width a. For 
such large values of n the Coulomb potential becomes 
a minor perturbation compared to to the infinite square 
well. By examining the diagonal matrix elements only 
for large values of n, one can derive 

(™a /a) 2 - 2^ ( 7 + In (2™)), (17) 

Eq a 

where 7 0.5772... is Euler's constant. This result is 
indicated with a curve in Fig. 2 and provides remark- 
ably accurate results (contrast with the ~ n 2 curve also 
shown), even for rather low values of n. While this ana- 
lytical result has nothing to do with the pure Coulomb 
potential it does provide an opportunity to illustrate first 
order perturbation theory 



TABLE II. Table 2. Results for a/a = 20. 





Ei/E 


E2/E0 


E 3 /E 


50 


-0.99428 


-0.24925 


-0.09951 


100 


-0.99920 


-0.24987 


-0.09979 


200 


-0.99989 


-0.24996 


-0.09983 


400 


-0.99999 


-0.24997 


-0.09984 



A summary of the effects of the square well cutoff and 
the matrix truncation size are best presented in tabular 
form, since the differences are so minute. Table 1 shows 
results as a function of the matrix size, n max , for a given 
a/ao = 10. In this instance the third eigenvalue is always 
positive, i.e. the square well is narrow enough that what 
would normally be the third bound state is pushed into 
the positive regime by the existence of the outer wall at 
r = a. Note that the first two bound states, tabulated in 
Table 1, do converge to a definite value as n max increases, 
but that this value is not necessarily the value pertain- 
ing to the Coulomb potential, without the embedding 
square well potential. In this case, this is especially true 
for E2/EQ, which should have a value of —0.25, but ac- 
tually converges to a value of —0.2256. If we didn't know 
beforehand that the expected values were Ei/Eq = — 1 
and E 2 /E = —1/4, then the way to check this is to 
increase the well width until we achieve convergence in 
the energies as a function of both a/ao and n max . Tables 
2 and 3 illustrate this process. It is clear that for suffi- 
ciently large a/ao the embedding potential plays no role 
(as is desirable) for a sufficiently large matrix size cutoff. 
In fact, for a specific cutoff, say n max = 50, then it is 
clear that as the size of the square well, a/ao, increases, 
the accuracy for a given energy level actually decreases. 
The reason for this is as stated earlier; for larger a/ao the 
same basis state (say, <f>5o(r)) has less spatial resolution 
than the 50 basis state for a smaller value of a/ao- 



TABLE III. Table 3. Results for a/a = 50. 



T^max 


E 1 /E 


E2/E0 


E3/E0 


50 


-0.94114 


-0.24207 


-0.10872 


100 


-0.98932 


-0.24864 


-0.11071 


200 


-0.99846 


-0.24981 


-0.11105 


400 


-0.99980 


-0.24997 


-0.11110 


800 


-0.99998 


-0.25000 


-0.11111 
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In summary we have illustrated how matrix mechan- 
ics, with a simple square well basis (i.e. a Fourier basis), 
can reproduce the bound state energies for the Coulomb 
potential. By extending the square well width, and ap- 
propriately increasing the matrix size, we converge to 
the known analytical results. We should also note that 
most software packages routinely return the eigenvec- 
tor along with the eigenvalue. For a given eigenvalue 
E n , this corresponds to a vector of coefficients Cm for 
m = l,2,3...n max , as in Eq. ([blj) . With these one can 
readily compute the corresponding radial wave function, 
R n (r) = u n (r)/r, where 



•-Jlaa(TEL) (18) 



and we have suppressed the index £ since we have fo- 
cussed on £ = for the Coulomb potential. 



square well disappears from the problem (as is desired) . 
We have also checked other states, including those with 
£ ^= 0, and found similar agreement. 



IV. FINITE SPHERICAL WELL 

The finite spherical well (see Fig. 1), with the simple 
potential, 



V sph {r) 



-Vq if < r < b, 
otherwise. 



(20) 



has the virtue that, at least for £ — 0, has matrix ele- 
ments that can be obtained analytically from Eq. (|10[) . 
They are 

H nm = S nm E° - — {g(n - m) - g(n + m)}, (21) 




FIG. 3. (color online) Numerical results (shown with sym- 
bols) for the radial wave functions (n = 1 and n — 2, both 
with I = 0) vs. the radial coordinate, r. Analytical results 
are also shown (with curves). Agreement is excellent. 



These are shown in Fig. 3 for the n = 1 and n 
states, along with the known analytical results, 



3/2 



3/2 



R w (r) = 2exp [-r/a ] 

i?2oM = -4(1 - ™ ) exp [-r/(2a<,)]; (19) 
V2 ^ &o 



the agreement of the numerical results is superb. Note 
that the scale is in units of the Bohr radius, and the 
decay is complete over a radial distance of about 10 x 
the Bohr radius. The right hand side of the embedding 
infinite square well is at a = 50ao, so this is well off- 
scale on this figure. Because both wave functions have 
decayed essentially to zero by this point, the embedding 



where as before E® — (Turn) 2 /[2moa 2 ], and the function 
g(n) is given by 



9{n) = 



sin (rnrb/a) 



(22) 



with the n = case simply given by l'Hopital's rule. 

Students can most readily work with this potential, 
and make comparisons with known analytical results. 
The analytical result in fact requires a graphical solu- 
tion of a transcendental equation^ 



tan (z + tt/2) 



/£0\2 

z ' 



1. 



(23) 



where z Q = n^/VojE 1 ^ and z = tt^(Vq + E)/E lb . Note 
that we use E° b = (hw) 2 j '(2mo& 2 ) as the energy scale, 
as the analytical solution does not utilize an embedding 
infinite potential well of width a. Solutions to Eq. ([23|) 
are easy to obtain if one simply solves for zq (and hence 
Vo) in terms of z (and hence E). Then, simple manipu- 
lation of Eq. (|23]l gives the result explicitly, from which 
a table can be readily constructed, and then E can be 
plotted vs. Vq. 

The spherical potential well can be used to illustrate 
the important principle that in three dimensions, a crit- 
ical depth Vqc is required to sustain a bound state, in 
contrast to the case in one or two dimensions. In fact, 
for a choice a — 106, the embedding potential will have 
no effect on the results, except as the potential depth is 
varied close to the critical potential. This is because as 
this occurs, the spatial extent of the bound state (for 
Vq ~ Voc) increases as Vq — > Vq c , until at some point, the 

embedding potential will 'aid' to push the bound state 
above zero energy for slightly larger Vq than would ac- 
tually occur without an embedding infinite square well. 

To discover this with our matrix approach, one would 
have to increase a/b, and determine the value of the po- 
tential for which a negative energy state no longer exists. 
A plot of these critical values of Vq, plotted versus b/a, 
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TABLE IV. Table 4. Ground state energies (A = 1) in units 
of -Bo, as compared with those of Ref. ((131 ). 



0.260 



.0' 



:] 




Is energy 


Ref. (13) 




0.10 


-0.814 116 


-0.814 116... 




0.20 


-0.653 617 


-0.653 617... 




0.40 


-0.396 752 


-0.396 752... 




0.60 


-0.212 271 


-0.212 272... 




0.80 


-0.089 408 


-0.089 409... 




1.00 


-0.020 552 


-0.020 572... 



0.255 



..0 



.0'' 



.•0 



0.250 



0.00 0.02 



0.04 0.06 

b/a 



0.10 



FIG. 4. (color online) A plot of Vbc, the critical value of Vb, 
below which a bound state no longer exists, vs. b/a, where b 
is the radius of the spherical potential, and a is the width of 
the embedding infinite square well potential. As a — > oo the 
influence of the embedding potential is eliminated; the blue 
line shows that the extrapolated value of Vbc = 1/4 agrees 
with the well known analytical result. 



is shown in Fig. 4. This clearly shows that as a — > oo, 
the critical value is Vo c /E® b — 0.25, in agreement with 
what is known analytically. We show this here to ad- 
dress a similar issue concerning the Yukawa potential in 
the next section. 



V. THE YUKAWA POTENTIAL 

We now examine an attractive Yukawa potential, 
which can be written as 



VVukM = -A 



-fir/a 



47re 



(24) 



where A allows us to adjust the strength of the interac- 
tion, and ao //j, is the screening length, written in units 
of the Bohr radius. Clearly, for A 1 and // — > 
we recover the Coulomb interaction discussed in Section 
III. The Yukawa potential has a shorter range than the 
Coulomb, and therefore has a finite number of bound 
states. In what follows we will focus exclusively on the 
£ = states, and thus have no need for the centrifugal 
term, although, as in the case of the unscreened Coulomb 
interaction, a study of £ ^ states poses no extra diffi- 
culty. 

Unlike the previous potentials discussed so far, there 
is no known analytical solution for the Yukawa poten- 
tial. Solutions either involve direct numerical solution 
of the Schrodinger differential equation,— or rely on so- 
phisticated numerical procedures centered around the 



variational method 9- ^ and perturbation theorj^; yet 
another numerical procedure uses an expansion tech- 
nique that uses special functions connected to Laguerre 
polynomials! 14 ' 15 In fact, our present methodology is 
similar in spirit to the so-called J-matrix method 16 
around which these Laguerre polynomial-based methods 
are developed. The virtue of the present approach is 
in its simplicity; the use of a simple numerical method 
based on straightforward matrix diagonalization in a ba- 
sis which consists of sine functions makes the work de- 
scribed below readily accessible to undergraduate stu- 
dents. 

Following the previous sections, we require the follow- 
ing matrix elements to be used in Eq. (1141) : 



= 5 nm (imao/a) 



- 2A—\K 1 (n + m)-K 1 (n-m)} 
a 

+£(£ + l)(J) 2 {L 2 (n + m) - L 2 (n ~ m)}, (25) 



with 

r , , , f 1 , [1 - cos (rmrx)]e-> iax / ao 

Ki(m) = / dx± ^ ^ , (26) 

Jo x 

where, as in the case for the Coulomb potential, we have 
used a unit of energy, Eq = ?i 2 /(2moag). Note that 
the important sections of numerical code needed to im- 
plement this calculation are included in an Appendix; 
setting fi = allows one to recover the results for the 
Coulomb interaction. 

The presence of an exponentially decaying factor in 
the K\ (m) integral does not make the (numerical) inte- 
gration any harder than in the Coulomb case; moreover, 
for most parameter choices converged results will be ob- 
tained without requiring a/a® to be excessively large. In 
Fig. 5 we show results for the energies of the s-states; 
symbols indicate previous results^ which are in excel- 
lent agreement with our own. Table 4 shows more digits 
for the ground state energies as a function of the screen- 
ing parameter fi, and indeed illustrates the remarkable 
accuracy of the results of Ref. (fl3l ). We have kept a/a 
and Umax fixed at 30 and 1800, respectively. This is the 
reason for the very slight deterioration in our values for 
larger /x; we can readily achieve further accuracy by in- 
creasing the infinite square well width, but at some point 
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FIG. 5. (color online) Energy levels for the Yukawa potential 
(with A = 1), as a function of the screening parameter, /u. 
The solid (red) , dotted (blue) and dashed (pink) curves show 
the Is, 2s, and 3s levels, respectively. Also shown (with 
symbols) are results from Ref. H3), with which we are in 
excellent agreement. Note the existence of critical values of 
He about which we will say more later. 
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FIG. 6. (color online) The ground state wave function vs. 
radial coordinate, r, for various values of the screening pa- 
rameter, /j,, while A has been increased to keep the ground 
state energy fixed, E\ — —Eq. As expected, the wave func- 
tion becomes increasingly more localized around the origin, 
as jj, increases, while keeping the binding energy constant. 



this would become prohibitively time-comsuming. As \x 
increases (for fixed A = 1), the bound state energy ap- 
proaches zero, and the wave function is more extended, 
and hence a large value of a is required to maintain the 
same accuracy. Note, however, that we are demonstrat- 
ing that we can achieve any desired accuracy; for exam- 
ple, with /i = 0.2, we require only n max ss 50 to achieve 
better than 0.1% accuracy in the energy. Figure 5 also 
illustrates graphically how quickly the infinite number of 
bound states become reduced to a very small number as 
p, increases. A critical value, fi c , beyond which no bound 
states exist, clearly exists near 1.2, about which we will 
say more below. 

As before one can obtain wave functions; when screen- 
ing is present, given the same energy (by increasing A 
as ji increases) the wave functions are more localized 
around the origin. Figure 6 shows a comparison of such 
wave functions; in each case we have adjusted the value 
of A to always maintain Ei(jj,)/Eq = — 1. The biggest 
impact occurs near the origin, as the screened wave func- 
tions have considerably more amplitude there. 

Conversely, we can examine how the ground state 
wave functions evolve with increasing fi with the strength 
maintained at A = 1. Then, the main effect will be that 
the energy approaches zero, so that the wave function 
is increasingly 'less bound' as /x increases. Thus, even 
though the interaction is more screened, and therefore 
shorter ranged, the wave function will spread out. Fig- 
ure 7 bears this out; as [i increases the wave functions be- 
come more extended, as expected. For \x w (x c = 1-1632 
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FIG. 7. (color online) The weighted ground state wave func- 
tion, ti(r) = rRio(r) vs. radial coordinate, r, for various val- 
ues of the screening parameter, jj,, while A is held constant 
at unity. Now as increases the binding energy decreases, 
so the wave function becomes more extended in space. Also 
shown is u(r) for fj, = fjb c ~ 1.1632 when a/ao = 50; now the 
wave function is on the verge of being delocalized over the 
entire space available, i.e. < r < a. 



(for the given value of a/ao = 50 used in Fig. 7) the 
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wave function extends over the entire space allowed by 
the infinite square well, and the shape is essentially that 
of a triangle in r: 

u(r) « B(l -e- Xr / ao )(l - r/a), (27) 

where B is determined through normalization. The first 
factor is required to ensure that u(r = 0) is zero. We 
find that A ps \[2[i c gives a remarkably good fit to the 
numerically attained wave function (it would be indis- 
tinguishable from the numerical result in Fig. 7). 
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FIG. 8. (color online) A plot of fi c , the critical value of fj,, 
above which a bound state no longer exists in the Yukawa 
potential, vs. oo/a, where ao is the Bohr radius, and a is 
the width of the embedding infinite square well potential. As 
in the case of the spherical well, as a —¥ oo the influence of 
the embedding potential is eliminated; the blue curve (fi c — 
1.1906 — 1.37ao/a) shows that the extrapolated value for the 
critical screening parameter is fi c ~ 1.1906, in agreement 
with previously established results. 

Finally, we wish to determine how to establish the 
critical value, fi c , above which no bound states exist. 
We proceed as in the case of the finite spherical well, 
and determine the critical value function of the 

infinite square well width, a/ag. The result is plotted 
in Fig. 8, and we determine, using the lowest two points 
(the most reliable) to extrapolate to a^/a — > 0, that /i c ps 
1.1906, in agreement with previous determinations ! 12 i 13 

VI. SUMMARY 

In this paper we have shown how one can use matrix 
mechanics, with the simplest of bases, to successfully 
obtain very accurate numerical results for the low-lying 
levels for essentially any three dimensional potential aris- 
ing from central forces that supports bound states. The 



mathematics required to do this is minimal, but stu- 
dents must be able to use any number of existing software 
packages to numerically diagonalize the resulting matrix. 
This skillset, non-existent a generation ago, is becoming 
increasingly useful for an undergraduate physics degree 
and beyond. 

Both bound state energies and wave functions can be 
readily obtained with the methodology described here. 
By increasing the size of the basis (and, if necessary, the 
size of the embedding infinite square well potential to 
accommodate the spatial spread of the bound state) one 
can achieve any desired accuracy. Results were demon- 
strated for the Coulomb, spherical well, and Yukawa po- 
tentials. We also addressed more difficult issues, such 
as the existence of critical parameters (well depth, or 
screening length) beyond which bound states cease to 
exist. We were able to reproduce the textbook result 
for the critical attractive potential for the spherical well, 
along with the not so well known result for the critical 
screening parameter in the case of the Yukawa potential. 

The properties of many other potentials used in the 
research literature are now accessible to undergraduates; 
studies of these potentials are suitable for assignments 
and/or projects. 
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Appendix A: A simple Fortran code to determine 
the eigenvalues and wave functions for the Yukawa 
Potential 

Most students will opt to use a high level language 
like MatLab, Mathematica, or Maple to solve problems 
as formulated in this paper. Only a few lines of code 
are required in this case. We have used C ++ and For- 
tran, and here we write down the key parts of the code 
required, in Fortran. We use a simple trapezoidal rule 
to evaluate the integrals required for the Hamiltonian 
matrix elements, and call upon two subroutines from 
Numerical Recipes^ to diagonalize the matrix. 

We use aa to designate the width of the well, while 
amu is the coefficient \x in the Yukawa potential, as writ- 
ten in Eq. (|2"^1) . We also use a coefficient amp to vary 
the strength of the Yukawa potential. If [i = then we 
have the Coulomb potential, and amp plays the role of 
Z, the nuclear charge. The key points of the code are as 
follows: 

c first get and save the needed integrals 
gg(0) = O.OdO 
gg2(0) = O.OdO 
do 11 n = l,2*nmax 
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gg(n) =0.0 
gg2(n) = 0.0 

do 226 iy = 2,nyy ! first term is taken 
care of separately below 

yy = yys + (iy - D*dyy 

term = (l.OdO - dcos (n*pi*yy) ) /yy 
gg(n) = gg(n) + term*dexp(-amu*aa*yy) 
gg2(n) = gg2(n) + term/yy 
226 continue 

c contribution from the origin 

gg(n) = dyy*(gg(n)+ O.OdO) ! nothing 
from the origin 

gg2(n) = dyy*(gg2(n) + 
. 5d0*0 . 5d0*n*n*pi*pi) ! 0.5 from trapezoidal 
rule 

11 continue 

The arrays gg(n) and gg2(n) contain the integrals 

Ki(n) and L2(n), respectively, as given in Eqs. (|2"B1) 

and (flU)) . The next bit of code constructs the matrix 

elements, h nmi as given in Eq. (fT5|) . 

c construct the needed matrix elements 
do 1 n = 1 , nmax 
do 2 m = l,n 

agl = 0.5d0*(gg(m+n) - gg(n-m)) 
ag2 = 0.5d0*(gg2(m+n) - gg2(n-m)) 
a(n,m) = 2 . 0d0*112*ag2/ (aa*aa) - 

4 . 0d0*amp*agl*zz/aa 

if (n.eq.m) a(n,m) = a(n,m) + 

(pi*n/aa) **2 

a(m,n) = a(n,m) 

2 continue 

1 continue 

The two dimensional array, a(n,m) now contains the 



Hamiltonian matrix given by Eq. (fT5)) . A call to the 
following two Numerical Recipes — routines then diago- 
nalizes and sorts the eigenvalues and eigenvectors, 
c these two routines, from Numerical Recipes, 
diagonalize the matrix a(n,m) 
call jacobi(a,nmax,dd,vv) 
call eigsrt(dd,vv,nmax) 
The eigenvalues are stored in the array dd(n) and the 
eigenvectors for the n th eigenvalue are stored in the two 
dimensional array, vv(m,n). That is, vv(m,n) is the 
coefficient for the m th basis state to the n th eigenvector. 
Thus, the next piece of code determines the ground state 
(n = 1) and the first excited state (n = 2) wave function 
as a function of x. 
c now for the wave function 
sq2 = dsqrt(2.0d0) 
do 5 ix = 1,5000 
xx = ix*0.0002d0 
bnsO = . OdO ! ground state 
bnsl = . OdO ! first excited state 
do 7 m = 1 , nmax 

bnsO = bnsO + vv(m, 1) *dsin(m*pi*xx) 
bnsl = bnsl + vv(m, 2) *dsin(m*pi*xx) 
7 continue 

bnsO = bns0*sq2 
bnsl = bnsl*sq2 
write(71,94)xx,bns0,bnsl 
5 continue 

94 f ormat(5x,f9.4,lx,f9.4,lx,f9.4) 

Similarly, the probabilities can be computed, and the 
total probability can be checked to sum to unity, as re- 
quired in the proper solution to the Schrodinger Equa- 
tion. 
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